export closure, small_generating_set

################################################################################
#
#  Find identity with respect to operation
#
################################################################################

# It is assumed that the elements have finite order with respect to op.
function find_identity(S, op)
  @assert length(S) > 0
  g = S[1]
  h = g
  while true
    hh = op(h, g)
    if hh == g
      return h
    end
    h = hh
  end
end

################################################################################
#
#  Small generating set
#
################################################################################

function small_generating_set(G, op)
  i = find_identity(G, op)
  return small_generating_set(G, op, i)
end

function small_generating_set(G, op, id)
  orderG = length(G)

  if length(G) == 1
    return G
  end

  function non_trivial_randelem()
    x = rand(G)
    while (x == id)
      x = rand(G)
    end
    return x
  end

  firsttry = 10
  secondtry = 20
  # First try one element
  for i in 1:firsttry
    gen = non_trivial_randelem()
    if length(closure([gen], op, id)) == orderG
      return [gen]
    end
  end

  for i in 1:secondtry
    gens = [non_trivial_randelem(), non_trivial_randelem()]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end

  # Now use that unconditionally log_2(|G|) elements generate G

  b = ceil(Int, log(2, orderG))
  @assert orderG <= 2^b

  j = 0
  while true
    if j > 2^20
      error("Something wrong with generator search")
    end
    j = j + 1
    gens = [non_trivial_randelem() for i in 1:b]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end
end

################################################################################
#
#  Computing closure under group operation
#
################################################################################

# It is assumed that S is nonempty and that the group generated by S under op
# is finite.
function closure(S, op)
  i = find_identity(S, op)
  return closure(S, op, i)
end

function closure(S, op, id)
  if length(S) == 0
    return [id]
  elseif length(S) == 1
    return _closing_under_one_generator(S[1], op, id)
  else
    return _closing_under_generators_dimino(S, op, id)
  end
end

function _closing_under_generators_naive(S, op, id)
  list = push!(copy(S), id)
  stable = false
  while !stable
    stable = true
    for g in list
      for s in S
        m = op(g, s)
        if !(m in list)
          push!(list, m)
          stable = false
        end
      end
    end
  end
  return list
end

function _closing_under_one_generator(x, op, id)
  elements = [x]
  y = x
  while !(y == id)
    y = op(y, x)
    push!(elements, y)
  end
  return elements
end

function _closing_under_generators_dimino(S, op, id)

  t = length(S)
  order = 1
  elements = [id]
  g = S[1]

  while g != id
    order = order +1
    push!(elements, g)
    g = op(g, S[1])
  end

  for i in 2:t
    if !(S[i] in elements)
      previous_order = order
      order = order + 1
      push!(elements, S[i])
      for j in 2:previous_order
        order = order + 1
        push!(elements, op(elements[j], S[i]))
      end

      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          s = S[k]
          elt = op(elements[rep_pos], s)
          if !(elt in elements)
            order = order + 1
            push!(elements, elt)
            for j in 2:previous_order
              order = order + 1
              push!(elements, op(elements[j], elt))
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  return elements
end

################################################################################
#
#  Multiplication table
#
################################################################################

# Construct multiplication table of G under op
function _multiplication_table(G, op)
  l = length(G)
  z = Matrix{Int}(l, l)
  for i in 1:l
    for j in 1:l
      p = op(G[i], G[j])
      for k in 1:l
        if p == G[k]
          z[i, j] = k
          break
        end
      end
    end
  end
  return z
end


################################################################################
#
#  Generic group given by multiplication table
#
################################################################################

# The elements are just 1..n and i * j = G.mult_table[i, j]
mutable struct GrpGen
  identity::Int
  order::Int
  mult_table::Array{Int, 2}
  gens::Array{Int, 1}

  function GrpGen(M::Array{Int, 2})
    z = new()
    z.mult_table = M
    z.identity = find_identity(z)
    return z
  end
end

struct GrpGenElem
  group::GrpGen
  i::Int
end

################################################################################
#
#  Compute generic group from anything
#
################################################################################

function generic_group(G, op)
  return GrpGen(_multiplication_table(G, op))
end

################################################################################
#
#  Parent
#
################################################################################

function parent(g::GrpGenElem)
  return g.group
end

################################################################################
#
#  Construct the ith element
#
################################################################################

function getindex(G::GrpGen, i::Int)
  return GrpGenElem(G, i)
end

function ==(g::GrpGenElem, h::GrpGenElem)
  return parent(g) == parent(h) && g.i == h.i
end

################################################################################
#
#  Order
#
################################################################################

function order(G::GrpGen)
  return size(G.mult_table, 1)
end

################################################################################
#
#  Order of an element
#
################################################################################

function order(g::GrpGenElem)
  k = 2
  h = g * g
  while h != g
    h = g * h
    k = k + 1
  end
  return k - 1
end

################################################################################
#
#  Identity
#
################################################################################

function find_identity(G::GrpGen)
  return _find_identity(G.mult_table)
end

function _find_identity(m::Array{Int, 2})
  return find_identity([1], (i, j) -> m[i, j])
end

function id(G::GrpGen)
  return GrpGenElem(G, G.identity)
end

################################################################################
#
#  Multiplication
#
################################################################################

function *(g::GrpGenElem, h::GrpGenElem)
  G = parent(g)
  return GrpGenElem(G, G.mult_table[g.i, h.i])
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(g::GrpGenElem)
  G = parent(g)
  m = G.mult_table
  i = g.i
  l = size(m, 1)
  ide = id(G)
  for j in 1:l
    if m[i, j] == G.identity
      return GrpGenElem(G, j)
    end
  end
  error("Not possible")
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, G::GrpGen)
  print(io, "Generic group with multiplication table\n")
  println(io, G.mult_table)
end

function Base.show(io::IO, g::GrpGenElem)
  print(io, "Element of generic group ($(g.i))\n")
end

################################################################################
#
#  Is abelian?
#
################################################################################

function isabelian(G::GrpGen)
  return defines_abelian_group(G.mult_table)
end

function defines_abelian_group(m)
  return issymmetric(m)
end

################################################################################
#
#  Isomorphism for specific groups
#
################################################################################

function isisomorphic_to_16T7(G::GrpGen)
  ordershape = [ (1, 1), (2, 3), (4, 12) ]
  # Use algorithm of Tarjahn (or something like that)
  # 16T7 has the following presentation:
  #
  #  $.1^4 = Id($)
  #  $.2^4 = Id($)
  #  $.3^2 = Id($)
  #  $.2^-1 * $.1^2 * $.2^-1 = Id($)
  #  $.1^-1 * $.2^-1 * $.1 * $.2^-1 = Id($)
  #  $.1^-1 * $.3 * $.1 * $.3 = Id($)
  #  $.2^-1 * $.3 * $.2 * $.3 = Id($)

  relations = [[(2, -1), (1, 2), (2, -1)],
               [(1, -1), (2, -1), (1, 1), (2, -1)], 
               [(1, -1), (3, 1), (1, 1), (3, 1)],
               [(2, -1), (3, 1), (2, 1), (3, 1)]]

  l = order(G)

  if l != 16
    return false
  end

  if isabelian(G)
    return false
  end

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  for i in 1:l
    g = G[i]
    o = order(g)
    if haskey(elements_by_orders, o)
      push!(elements_by_orders[o], g)
    else
      elements_by_orders[o] = [g]
    end
  end

  for (o, no) in ordershape
    if !haskey(elements_by_orders, o)
      return false
    else
      elts = elements_by_orders[o]
      if length(elts) != no
        return false
      end
    end
  end

  ide = id(G)

  img_of_gens = Iterators.product(elements_by_orders[4],
                                  elements_by_orders[4],
                                  elements_by_orders[2])

  for (g1, g2, g3) in img_of_gens
    g2inv = inv(g2)
    g1inv = inv(g1)
    #  $.2^-1 * $.1^2 * $.2^-1 = Id($)
    z = g2inv
    z = z * g1
    z = z * g1
    z = z * g2inv
    if z != ide
      continue
    end
    #  $.1^-1 * $.2^-1 * $.1 * $.2^-1 = Id($)
    z = g1inv
    z = z * g2inv
    z = z * g1
    z = z * g2inv
    if z != ide
      continue
    end
    #  $.1^-1 * $.3 * $.1 * $.3 = Id($)
    z = g1inv
    z = z * g3
    z = z * g1
    z = z * g3
    if z != ide
      continue
    end
    #  $.2^-1 * $.3 * $.2 * $.3 = Id($)
    z = g2inv
    z = z * g3
    z = z * g2
    z = z * g3
    if z != ide
      continue
    end

    cl = closure([g1, g2, g3], *, ide)

    if length(cl) < 16
      continue
    end
    return true
  end
  return false
end

function isisomorphic_to_8T5(G::GrpGen)
  l = order(G)

  if l != 8
    return false
  end

  # First check if it is abelian or not.
  if isabelian(G)
    return false
  end

  # Now G is D_8 (aka D_4 or 8T4) or Q_8
  # But D_8 has 5 elements of order 2, where as Q_8 has only 1 element of order 2.

  z = 0
  for i in 1:l
    g = G[i]
    o = order(g)
    if o == 2
      z += 1
    end
    if z > 1
      return false
    end
  end
  @assert z == 1
  return true
end
